#include "draft.h"
#include <stdlib.h>
#include "bss.h"
#include "mtp.h"

int max_hitsize;
static float best_score;        /* Best ordering score                */ 
static SeqCtg best_order[MAX_SEQCTG]; /* Ordered SeqCtgs with best 
                                       * ordering score               */ 
int batch = 0;

extern char bssFileLoad[];
extern BOOL next_mtp;
extern int num_sclone;
extern CLONE *seq_clone[];
extern char seq_name[];
extern int input_ctg, cur_contig, query_x, query_y;      
extern int max_seqctg, minolap;

extern void hit_CR();
extern int get_curctg(BSS_INST* pbsi);
extern void sort_seqctg_asce(SeqCtg **seqctgs, int size);
extern int read_seqctgs(SeqCtg *seq_ctgs);
extern void get_hits(HitList hits, SeqCtg *seq_ctg, int ctg, GArray *bes,BSS_INST* pbsi);
extern int in_bes(GArray *bes, BSS_HIT *t);
extern void sort_hitclone_start(HitClone **hit_clones, int size);

static void get_seqctginfo(SeqCtg *seq_ctgs, GArray *bes);
static void compute_scores(SeqCtg *seq_ctgs);
static void assign_recursive(GArray *bes, SeqCtg **s, int lev);
static void do_assign(GArray *bes, SeqCtg **s);
static void set_orientation(SeqCtg **s, GArray *bes, int included);
static void rc_recursive(GArray *bes, SeqCtg **s, int lev);
static void do_order(GArray *bes, SeqCtg **s);
static float order_medians(SeqCtg **s, GArray *bes); 
static void print_info(SeqCtg **s, GArray *bes);
static void print_order();
static float get_medians(int rc, HitClone **hits, int size);

static void sort_seqctg_median(SeqCtg** seq_ctgs, int size);
static void sort_seqctg_new(SeqCtg** seq_ctgs, int size);
static void sort_coord(int* coord, int size);

static gint seqctg_compare_median(const SeqCtg  **c1, const SeqCtg  **c2);
static gint seqctg_compare_new(const SeqCtg **c1, const SeqCtg **c2);
static gint coord_compare_asce(const int* c1, const int* c2);


GArray* hits_seqctg_ctg(BSS_INST* pbsi,int ctg, int sctg)
{
    GArray* ret;
    int i;
    BSS_HIT* hit;

    ret = g_array_new(FALSE,FALSE,sizeof(int));
    for (i = 0; i < arrayMax(pbsi->bss_hits); i++)
    {
        hit = arrp(pbsi->bss_hits,i,BSS_HIT);
        if (hit->data[BSS_CTG].i == ctg && hit->data[BSS_SEQCTG].i == sctg)
        {
            g_array_append_val(ret,i);
        } 
    }    
    return ret;
}

/*                        DEF : order_seqctgs
 *Order the sequence ctgs                                                */
void order_seqctgs() {
  BSS_INST bsi;
  gint ret;
  GArray **hit_clones;
  SeqCtg *seq_ctgs, **s;
  int i, j, hitsize;
  SeqCtg *seq_ctg;
  GArray *hit_list=NULL;
  HitClone *hit1, **best_hits[MAX_SEQCTG];
  GArray *bes;
  BESInfo *b1;

  hit_CR();

  if(strcmp(bssFileLoad, "") == 0) {
    printf("Please load BSS file first\n");
    return;
  }

  if (!read_bss_file(bssFileLoad,&bsi))
  {
    g_warning("Failed to  open %s", bssFileLoad);
    return;
  }

  next_mtp = FALSE;

  ret =  get_curctg(&bsi);
  if(ret == 2) {
    printf("The sequence is in multiple clones with differenct contigs.\n");
    for(i = 0; i < num_sclone; i++) 
      printf("Clone %s in contig %d\n", seq_clone[i]->clone, seq_clone[i]->ctg);
    printf("\n");
    return;
  } else if(ret == 1) {
    printf("There is no clone named %s\n\n", seq_name);
    return; 
  } else if(ret == 3) {
    printf("%s is in contig 0\n", seq_name);
    return;
  }

  if(input_ctg != -1 && input_ctg != cur_contig) {
    printf("Please select ctg%d.bss\n\n", cur_contig);
    return;
  }

  printf("\nQuery clone %s\n", seq_name);
  printf("Ctg%d", cur_contig);
  printf(" %d %d\n", query_x, query_y);

  seq_ctgs = (SeqCtg *)malloc(MAX_SEQCTG*sizeof(SeqCtg)); 
  if(seq_ctgs == NULL) 
    messcrash("ERROR---out of memory\n");

  max_seqctg = read_seqctgs(seq_ctgs);
 
  if(max_seqctg == 1) {
    printf("There is only one SeqCtg\n");
    return;
  } else if(max_seqctg == 0) {
    printf("There is no SeqCtgs in the BSS file.\n");
    return;
  }

  hit_clones = (GArray **)messalloc(sizeof(GArray *)*max_seqctg);
  bes = g_array_new(FALSE, FALSE, sizeof(BESInfo));
 
  max_hitsize = 0; 
  for(i = 0; i < max_seqctg; i++) { 
    seq_ctg = seq_ctgs+i; 
    hit_clones[i] = g_array_new(FALSE, FALSE, sizeof(HitClone));
    seq_ctg->hits = hit_clones[i];
    
    hit_list = hits_seqctg_ctg(&bsi,cur_contig,seq_ctg->seqctg); 
    if(hit_list != NULL) {
      get_hits(hit_list, seq_ctg, cur_contig, bes,&bsi);
      hitsize = seq_ctg->hits->len;
      if(hitsize > max_hitsize) max_hitsize = hitsize;
      g_array_free(hit_list,TRUE);
      hit_list = NULL;
    }
    for(j = 0; j < MAX_NEIGH; j++) {
      seq_ctg->dist_hits[j] = g_ptr_array_new();
      seq_ctg->neigh_hits[j] = g_ptr_array_new();
      seq_ctg->neigh_score[j] = 0;
      seq_ctg->dist_score[j] = 0;
      seq_ctg->neigh[j] = NULL;
      seq_ctg->distseqctg[j] = NULL;
    }
  }

  get_seqctginfo(seq_ctgs, bes);
  compute_scores(seq_ctgs);

  s = (SeqCtg **)malloc(max_seqctg*sizeof(SeqCtg *));
  if(s == NULL) 
    messcrash("ERROR---out of memory\n");

  for(i = 0; i < max_seqctg; i++) {
    s[i] = seq_ctgs+i;
  }

  sort_seqctg_asce(s, max_seqctg);
  print_info(s, bes);

  for(i = 0; i < max_seqctg; i++) {
    best_hits[i] = (HitClone **)malloc(max_hitsize*sizeof(HitClone *));
    if(best_hits[i] == NULL) 
      messcrash("ERROR---out of memory\n");
    best_order[i].hits_assigned = best_hits[i];
    best_order[i].count_assigned = 0;
  }

  best_score = 0.0;
  assign_recursive(bes, s, 0);
  print_order();

  /*Free memory allocated */
  for(i = 0; i < max_seqctg; i++) { 
    free(best_hits[i]);
    best_hits[i] = 0;
  }

  free(s);
  s = 0;

  for(i = 0; i < max_seqctg; i++) {
    for(j = 0; j < hit_clones[i]->len; j++) {
      hit1 = &g_array_index(hit_clones[i], HitClone, j);
      g_ptr_array_free(hit1->neigh_hits, TRUE); 
      g_ptr_array_free(hit1->dist_hits, TRUE); 
    }
    g_array_free(hit_clones[i], TRUE);
  }
  messfree(hit_clones);

  for(i = 0; i < max_seqctg; i++) { 
    seq_ctg = seq_ctgs+i; 
    for(j = 0; j < MAX_NEIGH; j++) {
       g_ptr_array_free(seq_ctg->dist_hits[j], TRUE);
       g_ptr_array_free(seq_ctg->neigh_hits[j], TRUE);
    }
  }
 
  for(i = 0; i < bes->len; i++) {
    b1 = &g_array_index(bes, BESInfo, i);
    g_ptr_array_free(b1->seqctg_list, TRUE);
  } 
  g_array_free(bes, TRUE);
  free(seq_ctgs);
  seq_ctgs = 0;
}

/*                        DEF : assign_recursive
 * Recursively assign all BES uniquely to a SeqCtg */
static void assign_recursive(GArray *bes, SeqCtg **s, int lev) {
  GPtrArray *list;
  BESInfo *b;
  int i, len;
  SeqCtg *seq_ctg;

  if(lev == bes->len) {
    do_assign(bes, s);
    return;
  }

  b = &g_array_index(bes, BESInfo, lev);
  list = b->seqctg_list;
  if(b->in_neigh) len = 1;
  else len = list->len;

  for(i = 0; i < len; i++) {
    seq_ctg = g_ptr_array_index(list, i);
    b->seqctg_assigned = seq_ctg->seqctg;
    assign_recursive(bes, s, lev+1);
  } 
}

/*                        DEF : do_assign 
 * Compute the orientation, median, and related info of 
 * SeqCtgs based on the current assignment             */ 
static void do_assign(GArray *bes, SeqCtg **s) {
  HitClone **hits[MAX_SEQCTG];
  HitClone *h;
  BSS_HIT *t;
  float xend_rc_y, xend_rc_n, yend_rc_y, yend_rc_n, score;
  BESInfo *b;
  int index;
  int i, j;
  SeqCtg *seq_ctg;

  for(i = 0; i < max_seqctg; i++) { 
    seq_ctg = s[i]; 
    seq_ctg->included = 0;
    seq_ctg->median = -1;
    seq_ctg->oret = UNKNOWN;
    seq_ctg->oret_score = -1;
    seq_ctg->hits_assigned = NULL;
    seq_ctg->count_assigned = 0;
    seq_ctg->rc_n_score = 0;
    seq_ctg->rc_y_score = 0;

    hits[i] = (HitClone **)malloc(seq_ctg->hits->len*sizeof(HitClone *));
    if(hits[i] == NULL) 
      messcrash("ERROR---out of memory\n");

    xend_rc_y = xend_rc_n = yend_rc_y = yend_rc_n = 0.0;

    for(j = 0; j < seq_ctg->hits->len; j++) { 
      h = &g_array_index(seq_ctg->hits, HitClone, j);
      t = h->hit;

      index = in_bes(bes, t);
      b = &g_array_index(bes, BESInfo, index);

      if(seq_ctg->seqctg == b->seqctg_assigned || b->in_neigh) {
        hits[i][seq_ctg->count_assigned] = h;
        score = 1.0/h->count;
        seq_ctg->count_assigned++;
        if(h->end_in == XIN) {
          if(bss_is_rc(t) == 1) 
            xend_rc_y += score;
          else
            xend_rc_n += score; 
        }  
        else if(h->end_in == YIN) { 
          if(bss_is_rc(t) == 1) 
            yend_rc_y += score;
          else
            yend_rc_n += score; 
        }
      }
    }
    seq_ctg->rc_n_score = xend_rc_n + yend_rc_y;
    seq_ctg->rc_y_score = xend_rc_y + yend_rc_n;

    if(seq_ctg->rc_y_score > 0.0 || seq_ctg->rc_n_score > 0.0) { 
      /*Find SeqCtgs that can be ordered */ 
      seq_ctg->included = 1;
      if(seq_ctg->rc_y_score == 0.0) { 
        seq_ctg->oret = RCN;
        seq_ctg->oret_score = seq_ctg->rc_n_score;
      }  else if(seq_ctg->rc_n_score == 0.0) {
        seq_ctg->oret = RCY;
        seq_ctg->oret_score = seq_ctg->rc_y_score;
      } 
    } else {
      if(seq_ctg->neigh_count > 0 || seq_ctg->dist_count > 0) 
        seq_ctg->included = 3;
    } 
    seq_ctg->hits_assigned = hits[i];
  }

  /*Set the orienations of neighbors or distant seqctgs of the seqctgs 
   *whose orientaion can be determined */ 
  set_orientation(s, bes, 1);

  /*Set the orienations of neighbors or distant seqctgs of the seqctgs 
   *whose orientaion is set by its other neighbors or distant seqctgs */ 
  set_orientation(s, bes, 2);

  rc_recursive(bes, s, 0);

  for(i = 0; i < max_seqctg; i++) { 
    free(hits[i]);
    hits[i] = 0;
  }
}

static void set_orientation(SeqCtg **s, GArray *bes, int included) {
  int i, j, k, l, seqctg1, seqctg2;
  BESInfo *b1, *b2;
  BOOL found = FALSE; 
  SeqCtg *seq_ctg, *neigh_seqctg, *dist_seqctg;
  int neigh_count, dist_count;
  HitClone *hit1, *hit2;
  BSS_HIT *t1, *t2;

  for(i = 0; i < max_seqctg; i++) {
    seq_ctg = s[i]; 
    seqctg1 = seq_ctg->seqctg;
    neigh_count = seq_ctg->neigh_count;
    dist_count = seq_ctg->dist_count;

    found = FALSE;
    if(seq_ctg->included == included) {
      if(neigh_count > 0) {
        for(j = 0; j < neigh_count; j++) {
          neigh_seqctg = seq_ctg->neigh[j];
          seqctg2 = neigh_seqctg->seqctg;

          if(neigh_seqctg->oret == UNKNOWN && 
            (neigh_seqctg->included == 0 || neigh_seqctg->included == 3)) {
            for(k = 0; k < seq_ctg->neigh_hits[j]->len; k++) {
              hit1 = g_ptr_array_index(seq_ctg->neigh_hits[j], k);
              hit2 = g_ptr_array_index(neigh_seqctg->neigh_hits[j], k);
              t1 = hit1->hit;
              t2 = hit2->hit;

              if(bss_is_rc(t1) != bss_is_rc(t2)) 
                neigh_seqctg->oret = seq_ctg->oret == RCN ? RCY : RCN; 
              else 
                neigh_seqctg->oret = seq_ctg->oret;
              neigh_seqctg->oret_score = seq_ctg->oret_score;
              found = TRUE;
              neigh_seqctg->included = 2;
              break;
            }
          } 
        } 
      } else if(dist_count > 0) {
        for(j = 0; j < dist_count; j++) {
          dist_seqctg = seq_ctg->distseqctg[j];
          seqctg2 = dist_seqctg->seqctg;

          if(dist_seqctg->oret == UNKNOWN &&  
            (dist_seqctg->included == 0 || dist_seqctg->included == 3)) {
            for(k = 0; k < seq_ctg->dist_hits[j]->len; k++) {
              hit1 = g_ptr_array_index(seq_ctg->dist_hits[j], k);
              b1 = &g_array_index(bes, BESInfo, hit1->bes_index);
              if(b1->seqctg_assigned == seqctg1) { 
                for(l = 0; l < hit1->dist_hits->len; l++) {
                  hit2 = g_ptr_array_index(hit1->dist_hits, l);
                  if(hit2->hit->data[BSS_SEQCTG].i == seqctg2) {
                    b2 = &g_array_index(bes, BESInfo, hit2->bes_index);

                    if(b2->seqctg_assigned == seqctg2) {  
                      t1 = hit1->hit;
                      t2 = hit2->hit;

                      if(bss_is_rc(t1) != bss_is_rc(t2))
                        dist_seqctg->oret = seq_ctg->oret;
                      else 
                        dist_seqctg->oret = seq_ctg->oret == RCN ? RCY : RCN;

                      dist_seqctg->oret_score = seq_ctg->oret_score;
                      found = TRUE;
                      dist_seqctg->included = 2;
                      break;
                    }
                  }
                }
                if(found) break;
              }
            }
          }
        }
      }
    }
  }
}

/*                        DEF : rc_recursive 
 * Recursively set the orienation of the SeqCtg if its  
 * orientation can not be determined                    */ 
static void rc_recursive(GArray *bes, SeqCtg **s, int lev) {
  orientation rc[2];
  float rc_score[2];
  int i, len, flag = 0;

  if(lev == max_seqctg) {
    do_order(bes, s);
    return;
  }

  rc[0] = RCN;
  rc[1] = RCY;

  if((s[lev]->rc_n_score > 0.0 && s[lev]->rc_y_score > 0.0) || s[lev]->included == 3) {
    flag = 1;
    len = 2;
    if(s[lev]->included != 3) {
      rc_score[0] = s[lev]->rc_n_score;
      rc_score[1] = s[lev]->rc_y_score;
    } else {
      rc_score[0] = 0;
      rc_score[1] = 0;
    }
  } else 
    len = 1; 

  for(i = 0; i < len; i++) {
    if(flag) {
      s[lev]->oret = rc[i];
      s[lev]->oret_score = rc_score[i];
    }
    rc_recursive(bes, s, lev+1);
  }
}

/*                        DEF : do_order 
 * Order the SeqCtgs based on current assignment and its 
 * orientation                                           */ 
static void do_order(GArray *bes, SeqCtg **s) { 
  float total_score;
  HitClone **hits;
  int i, j, count;
  SeqCtg *seq_ctg;

  for(i = 0; i < max_seqctg; i++) {
    seq_ctg = s[i]; 
    /*compute the median of clone coordinates */
    if(seq_ctg->included) {
      hits = seq_ctg->hits_assigned;
      count = seq_ctg->count_assigned; 
      if(seq_ctg->oret == RCN) 
        seq_ctg->median = get_medians(0, hits, count);
      else if(seq_ctg->oret == RCY) 
        seq_ctg->median = get_medians(1, hits, count);
    }
  }

  total_score = order_medians(s, bes);
  if(total_score > best_score) {
    best_score = total_score;
   
    for(i = 0; i < max_seqctg; i++) {
      best_order[i].seqctg = s[i]->seqctg;
      best_order[i].new_seqctg = s[i]->new_seqctg;
      best_order[i].median = s[i]->median;
      best_order[i].oret = s[i]->oret;
      best_order[i].oret_score = s[i]->oret_score;
      best_order[i].included = s[i]->included;
      best_order[i].count_assigned = 0;
      hits = best_order[i].hits_assigned; 
      for(j = 0; j < s[i]->count_assigned; j++) {
        hits[j] = s[i]->hits_assigned[j];
        hits[j]->best_in = s[i]->hits_assigned[j]->contained_in;
      }
      best_order[i].count_assigned = j;
    }
  }
}

/*                         DEF: get_seqctginfo
 *Gather all information about the hits of SeqCtgs */
static void get_seqctginfo(SeqCtg *seq_ctgs, GArray *bes) {
  SeqCtg *s1, *s2;
  GArray *hits1, *hits2;
  HitClone *h1, *h2;
  BSS_HIT *t1, *t2;
  int i, j, k, l, m, count = 0, index;  
  int index1, index2, hit_len1, hit_len2;
  int ncount1, ncount2, dcount1;
  int seqctg1, seqctg2;
  HitClone *tmp_hits[MAX_SEQCTG];
  CLONE *clp1, *clp2;
  BESInfo *b;
 
  for(i = 0; i < max_seqctg; i++) {
    s1 = seq_ctgs+i;
    hits1 = s1->hits;
    ncount1 = s1->neigh_count;
    dcount1 = s1->dist_count;
    seqctg1 = s1->seqctg;

    /*Traverse all hits in SeqCtg1 */
    for(j = 0; j < hits1->len; j++) {
      h1 = &g_array_index(hits1, HitClone, j);
      t1 = h1->hit; 
      clp1 = h1->clp;      

      /*check all hits in other seqctgs */ 
      for(k = 0; k < max_seqctg; k++) {
        if(k == i) continue;
        s2 = seq_ctgs+k; 
        hits2 = s2->hits;
        ncount2 = s2->neigh_count;
        seqctg2 = s2->seqctg;

        /*Traverse all hits in SeqCtg2 */
        for(l = 0; l < hits2->len; l++) {
          h2 = &g_array_index(hits2, HitClone, l);
          t2 = h2->hit;
          clp2 = h2->clp;

          /*if bes1 and bes2 hit the same clone */
          if(!strcmp(clp1->clone, clp2->clone)) { 
         
            /*if bes1 and bes2 are identical */
            if(!strcmp(t1->data[BSS_BES].str, t2->data[BSS_BES].str)) {
              /* if the bes is not checked       */
              if(!h2->checked && !h1->checked) {

                /* If the hit is at the end of both SeqCtgs and the sum of hit 
                 * length is less than the BES length, the two SeqCtgs are neighbors*/
                if(h1->end != 0 && h2->end != 0) { 
                  hit_len1 = t1->data[BSS_QEND].i - t1->data[BSS_QSTART].i + 1;
                  hit_len2 = t2->data[BSS_QEND].i - t2->data[BSS_QSTART].i + 1;
                  if((hit_len1 + hit_len2) <= t1->data[BSS_TLEN].i) {
                    index1 = -1;
                    index2 = -1;
                    for(m = 0; m < ncount1; m++) {
                      if(seqctg2 == s1->neigh[m]->seqctg) 
                        index1 = m;
                    }

                    if(index1 < 0) {  
                      s1->neigh[ncount1] = s2;
                      g_ptr_array_add(s1->neigh_hits[ncount1], h1);
                      s1->neigh_count++;
                      ncount1++;
                      index = in_bes(bes, h1->hit);
                      b = &g_array_index(bes, BESInfo, index);
                      b->in_neigh = 1;
                    } else 
                      g_ptr_array_add(s1->neigh_hits[index1], h1);
                    g_ptr_array_add(h1->neigh_hits, h2);

                    for(m = 0; m < ncount2; m++) {
                      if(seqctg1 == s2->neigh[m]->seqctg) 
                        index2 = m;
                    }
                    
                    if(index2 < 0) {  
                      s2->neigh[ncount2] = s1;
                      g_ptr_array_add(s2->neigh_hits[ncount2], h2);
                      s2->neigh_count++;
                      ncount2++;
                    } else 
                      g_ptr_array_add(s2->neigh_hits[index2], h2);
                    g_ptr_array_add(h2->neigh_hits, h1);
                  }
                } 
                h1->count++;
                tmp_hits[count++] = h2;
                h2->checked = 1;
              } 
            } else { /*If bes1 and bes2 are two ends of the same clone */
              index1 = -1;
              for(m = 0; m < dcount1; m++) {
                if(seqctg2 == s1->distseqctg[m]->seqctg) 
                  index1 = m;
              }
              if(index1 < 0) {  
                s1->distseqctg[dcount1] = s2;
                g_ptr_array_add(s1->dist_hits[dcount1], h1);
                s1->dist_count++;
                dcount1++;
              } else  
                g_ptr_array_add(s1->dist_hits[index1], h1);
              g_ptr_array_add(h1->dist_hits, h2);
            } 
          } 
        }  /*end of for of l */
      } /*end of for of k */   

      for(k = 0; k < count; k++) 
        tmp_hits[k]->count = h1->count;

      count = 0;
      h1->checked = 1;
    } /*end of for of j */
  } /* end of for of i */
} 

/*                         DEF: compute_scores 
 *Compute the scores                             */
static void compute_scores(SeqCtg *seq_ctgs) {
  int start, end, i, j, k;
  SeqCtg *s;
  GArray *hits;
  int ncount, dcount;
  HitClone *h;
  BSS_HIT *t;
  CLONE *clp;
  float score;

  start = query_x-minolap;
  end = query_y+minolap;

  for(i = 0; i < max_seqctg; i++) {
    s = seq_ctgs+i; 
    hits = s->hits;
    ncount = s->neigh_count;
    dcount = s->dist_count;
   
    /*Compute neighbor scores */ 
    if(ncount <= 2) {
      for(j = 0; j < ncount; j++) 
        s->neigh_score[j] = s->neigh_hits[j]->len;
    } else {
      for(j = 0; j < ncount; j++) 
        s->neigh_score[j] = s->neigh_hits[j]->len/(float)ncount; 
    }  
     
    /*Compute distant seqctg scores */ 
    for(j = 0; j < dcount; j++) {
      for(k = 0; k < s->dist_hits[j]->len; k++) {
        h = g_ptr_array_index(s->dist_hits[j], k);
        s->dist_score[j] += 1.0/h->count;
      }
    }
    
    for(j = 0; j < hits->len; j++) {
      h = &g_array_index(hits, HitClone, j);
      t = h->hit; 
      clp = h->clp;
      score = 1.0/h->count;
      
      if(h->parent_hit) 
        s->parent_score = 1.0/h->count;
    }
  }
}                
         
/*                        DEF : order_medians
 *Order the sequence ctgs based on its coordinate median           */
static float order_medians(SeqCtg **s, GArray *bes) {
  int i, j, k, l, count = 0, ncount, dcount, max;
  BSS_HIT *t1, *t2;
  SeqCtg *seq_ctg, *dist_seqctg, *neigh_seqctg;
  HitClone *hit1, *hit2;
  float total_score=0.0;
  int seqctg1, seqctg2;
  BESInfo *b1, *b2;

  sort_seqctg_median(s, max_seqctg);
  count = 1;
  for(i = 0; i < max_seqctg; i++) {
    seq_ctg = s[i]; 
    if(seq_ctg->median == -1)  {
      seq_ctg->oret = UNKNOWN;
      seq_ctg->included = 0;
      seq_ctg->new_seqctg = max_seqctg;
    }
    if(seq_ctg->included) 
      seq_ctg->new_seqctg = count++;
  }
  sort_seqctg_new(s, max_seqctg);

  max = count -1;
  for(i = 0; i < max_seqctg; i++) {
    seq_ctg = s[i]; 
    seqctg1 = seq_ctg->seqctg;
    ncount = seq_ctg->neigh_count;
    dcount = seq_ctg->dist_count;
    
    if(seq_ctg->included == 1 || seq_ctg->included == 3) {
      /* add the orienation score to the total ordering score */
      total_score += seq_ctg->oret_score;

      for(j = 0; j< ncount; j++) { 
        neigh_seqctg = seq_ctg->neigh[j]; 
        seqctg2 = neigh_seqctg->seqctg;

        /* If two possible neighbors or distant SeqCtgs are actually 
         * neighbors or distant SeqCtgs in the order, add their scores to 
         * the total score */
        if(neigh_seqctg->included == 1 || neigh_seqctg->included == 3) {
          for(k = 0; k < seq_ctg->neigh_hits[j]->len; k++) {
            hit1 = g_ptr_array_index(seq_ctg->neigh_hits[j], k);
            hit2 = g_ptr_array_index(neigh_seqctg->neigh_hits[j], k);
            t1 = hit1->hit;
            t2 = hit2->hit;

            if(abs(neigh_seqctg->new_seqctg - seq_ctg->new_seqctg) == 1) { 
              if((seq_ctg->oret == neigh_seqctg->oret && bss_is_rc(t1) == bss_is_rc(t2)) ||
                 (seq_ctg->oret != neigh_seqctg->oret && bss_is_rc(t1) != bss_is_rc(t2)))  {
                total_score += seq_ctg->neigh_score[j];
              }
            }
          }
        } else if(neigh_seqctg->included == 2) {
          if(abs(neigh_seqctg->new_seqctg - seq_ctg->new_seqctg) == 1) 
            total_score += seq_ctg->neigh_score[j]; 
        }
      }

      for(j = 0; j< dcount; j++) { 
        dist_seqctg = seq_ctg->distseqctg[j];
        seqctg2 = dist_seqctg->seqctg;

        if(dist_seqctg->included == 1 || dist_seqctg->included == 3) {
          for(k = 0; k < seq_ctg->dist_hits[j]->len; k++) {
            hit1 = g_ptr_array_index(seq_ctg->dist_hits[j], k);
            b1 = &g_array_index(bes, BESInfo, hit1->bes_index);

            if(b1->seqctg_assigned == seqctg1) {
              for(l = 0; l < hit1->dist_hits->len; l++) {
                hit2 = g_ptr_array_index(hit1->dist_hits, l);
                if(hit2->hit->data[BSS_SEQCTG].i == seqctg2) {
                  b2 = &g_array_index(bes, BESInfo, hit2->bes_index);

                  if(b2->seqctg_assigned == seqctg2) {  
                    t1 = hit1->hit;
                    t2 = hit2->hit;
                    if(abs(dist_seqctg->new_seqctg - seq_ctg->new_seqctg) > 1)  {
                      if((seq_ctg->oret == dist_seqctg->oret && bss_is_rc(t1) != bss_is_rc(t2)) ||
                         (seq_ctg->oret != dist_seqctg->oret && bss_is_rc(t1) == bss_is_rc(t2)))  {
                        total_score += seq_ctg->dist_score[j];
                      }
                    }
                  }
                }
              }
            }
          }
        } else if (dist_seqctg->included == 2) {
          if(abs(dist_seqctg->new_seqctg - seq_ctg->new_seqctg) > 1)  
            total_score += seq_ctg->dist_score[j];
        }
      }
    }

    /* If the possible end SeqCtgs are really at the either end, add the score to 
     * the total ordering score                                                 */
    if(seq_ctg->parent_score > 0) 
      if(seq_ctg->new_seqctg == 1 || seq_ctg->new_seqctg == max)
        total_score += seq_ctg->parent_score;
  }

  return total_score;
}

/*                        DEF : print_info 
 * Print all info about the SeqCtgs                */
static void print_info(SeqCtg **s, GArray *bes) {
  int i, j, k, count=0, ncount, dcount;
  char end;
  HitClone **hits;
  char str_format[100];
  char tmp_rc;
  BSS_HIT *t;
  HitClone *h;
  SeqCtg *seq_ctg, *tmp_seqctg;
  BESInfo *b;

  for(i = 0; i < max_seqctg; i++) {
    seq_ctg = s[i];
    hits = (HitClone **)malloc(seq_ctg->hits->len*sizeof(HitClone *));
    if(hits == NULL) 
      messcrash("ERROR---out of memory\n");

    count = 0;
    for(j = 0; j < seq_ctg->hits->len; j++) 
      hits[count++] = &g_array_index(seq_ctg->hits, HitClone, j);
    sort_hitclone_start(hits, count);

    if(count > 0) {
      printf("Seqctg %d:\n", seq_ctg->seqctg);
      if(seq_ctg->parent_score > 0)  
        printf("Contain parent BES; score: %.2f\n", seq_ctg->parent_score);
    }
    if(seq_ctg->hits->len != 0)  
      printf("BES               RC Start      End        x     y   HitLen%% Identity #Hits SeqctgList\n");
    
    for(j = 0; j < count; j++) {
      h = hits[j];
      t = h->hit;
      tmp_rc = bss_is_rc(t) ? 'y' : 'n';
      if(h->end) end = '*';
      else end = ' ';
    
      if(h->end_in == XIN)  
         strcpy(str_format, "%-18s%-3c%-10u%-10u -%-5u %-4u%5d%%%c %5d%%%5u    ");
      else if (h->end_in == YIN) 
         strcpy(str_format, "%-18s%-3c%-10u%-10u  %-5u-%-4u%5d%%%c %5d%%%5u    ");
      else 
         strcpy(str_format, "%-18s%-3c%-10u%-10u -%-5u-%-4u%5d%%%c %5d%%%5u    ");
      printf(str_format, 
            t->data[BSS_BES].str, tmp_rc, t->data[BSS_QSTART].i, t->data[BSS_QEND].i, 
            h->clp->x,  h->clp->y, (int)h->pct_hitlen, end, t->data[BSS_IDENTITY].i, h->count);

      b = &g_array_index(bes, BESInfo, h->bes_index);
      for(k = 0; k < b->seqctg_list->len; k++) {
        tmp_seqctg = g_ptr_array_index(b->seqctg_list, k);
        if(k == 0) 
          printf("Seqctg%d", tmp_seqctg->seqctg);
        else 
          printf(", %d", tmp_seqctg->seqctg);
      }
      printf("\n"); 
    }
    if(count != 0)  
      printf("\n");

    ncount = seq_ctg->neigh_count;
    dcount = seq_ctg->dist_count;

    if(ncount > 0) {
      printf("%d Neighbors: \n", ncount);
      printf("Neighbor     Score  #Hits  HitList\n");
      for(j = 0; j< ncount; j++) { 
        printf("Seqctg%2d     %.2f%6d    ", seq_ctg->neigh[j]->seqctg, 
              seq_ctg->neigh_score[j], seq_ctg->neigh_hits[j]->len);
        for(k = 0; k < seq_ctg->neigh_hits[j]->len; k++) {
          h = g_ptr_array_index(seq_ctg->neigh_hits[j], k);
          t = h->hit;
          if(k == 0) 
            printf("%s", t->data[BSS_BES].str);
          else 
            printf(", %s", t->data[BSS_BES].str);
        }
        printf("\n");
      }
      printf("\n");
    }

    if(dcount > 0) {
      printf("%d Distant Seqctgs: \n", dcount);
      printf("DistSeqctg   Score  #Hits  HitList\n");
      for(j = 0; j< dcount; j++) { 
        printf("Seqctg%2d     %.2f%6d    ", seq_ctg->distseqctg[j]->seqctg, 
              seq_ctg->dist_score[j], seq_ctg->dist_hits[j]->len);
        for(k = 0; k < seq_ctg->dist_hits[j]->len; k++) {
          h = g_ptr_array_index(seq_ctg->dist_hits[j], k);
          t = h->hit;
          if(k == 0) 
            printf("%s", t->data[BSS_BES].str);
          else 
            printf(", %s", t->data[BSS_BES].str);
        }
        printf("\n");
      }
      printf("\n");
    }

    free(hits);
    hits=0;
    printf("\n");
  }
}

/*                         DEF: print_order  
 *Print the ordered SeqCtgs                   */
static void print_order() {
  int i, j, count=0;
  FILE *fp;
  char end;
  char str_tmp[10], str_output[100], str_format[100];
  char tmp_rc;
  BSS_HIT *t;
  HitClone *h;
  SeqCtg *seq_ctg;
  HitClone **hits;

  strcpy(str_output, " (number of hits: ");
  if(batch) {
    fp = fopen("seqctg_list", "a");
    if(fp == NULL) 
      fp = fopen("seqctg_list", "w");
    fprintf(fp, ">%s.seq\n", seq_name);
    fprintf(fp, "SegCtg");
  } 
  printf("Best order score %.2f\n", best_score);

  for(i = 0; i < max_seqctg; i++) {
    seq_ctg = best_order+i;
    hits = seq_ctg->hits_assigned;
    count = seq_ctg->count_assigned;
    sort_hitclone_start(hits, count);
    if(seq_ctg->median != -1)  
      printf("\nSeqCtg %d Median: %.2f\n", seq_ctg->seqctg, seq_ctg->median); 
    else if(seq_ctg->included) 
      printf("\nSeqCtg %d\n", seq_ctg->seqctg); 
    
    if(seq_ctg->included) {
      printf("Orientation: ");
      if(seq_ctg->oret == RCN)  
        printf("right orienation; score: %.2f\n", seq_ctg->oret_score);
      else  
        printf("reverse complemented; score: %.2f\n", seq_ctg->oret_score);
    }

    if(batch) { 
      if(seq_ctg->included) {
        fprintf(fp, " %d", seq_ctg->seqctg);
        sprintf(str_tmp, "%d ", count);
        strcat(str_output, str_tmp);
      }
    } 
    if(seq_ctg->included) {
      if(count != 0)  
        printf("BES              RC Start      End        x        y      HitLen%%    Identity\n");
      
      for(j = 0; j < count; j++) {
        h = hits[j];
        t = h->hit;
        tmp_rc = bss_is_rc(t) ? 'y' : 'n';
        if(h->end) end = '*';
        else end = ' ';

        if(h->best_in == XIN)  
          strcpy(str_format, "%-18s%-3c%-10u%-10u -%-8u %-8u%3d%%%c     %5u%%\n");
        else 
          strcpy(str_format, "%-18s%-3c%-10u%-10u  %-8u-%-8u%3d%%%c     %5u%%\n");
        printf(str_format, t->data[BSS_BES].str, tmp_rc, t->data[BSS_QSTART].i, 
              t->data[BSS_QEND].i, h->clp->x, h->clp->y, (int)h->pct_hitlen, end, 
                t->data[BSS_IDENTITY].i);
      }
      if(count != 0)  
        printf("\n");
    }
  }
  
  if(batch) {
    strcat(str_output, ")");
    fprintf(fp, "%s\n", str_output); 
    fclose(fp);
  }
}

/*                         DEF: get_medians  
 *Get the median of relevent coordinates for a SeqCtg      */
static float get_medians(int rc, HitClone **hits, int size) {
    int *coordinates;
    int i, count=0;
    float median;
    HitClone *h;
    BSS_HIT *t;
    CLONE *clp;

    coordinates = (int *)malloc(size*sizeof(int));
    if(coordinates == NULL) 
      messcrash("ERROR---out of memory\n");

    for(i = 0; i < size; i++) {
      h = hits[i];
      t = h->hit; 
      clp = h->clp; 
     
      if(h->end_in == XIN) { 
        coordinates[count++] = clp->x;
        h->contained_in = XIN;
      } else if (h->end_in == YIN) {
        coordinates[count++] = clp->y; 
        h->contained_in = YIN;
      } else { 
        if(bss_is_rc(t) == rc)  { 
          coordinates[count++] = clp->x;
          h->contained_in = XIN;
        }
        else {
          coordinates[count++] = clp->y;
          h->contained_in = YIN;
        }
      }
    }
    
    sort_coord(coordinates, count);
    if(count != 0) {
      if(count%2 == 0)  
        median = (float)(coordinates[count/2-1]+coordinates[count/2])/2;
      else  
        median = coordinates[count/2];
    } else 
        median = -1;
    
    free(coordinates);
    coordinates = 0;
    return median;
}

/*                   DEF: sort_seqctg_median 
 *Sort the sequence contigs in ascending order of medians*/
static void sort_seqctg_median(SeqCtg** seq_ctgs, int size)
{
  qsort(seq_ctgs, size, sizeof(SeqCtg *),
    (int (*)(const void*, const void*))seqctg_compare_median);
}

/*                   DEF: sort_seqctg_new 
 *Sort the sequence contigs in ascending order of ordered SeqCtg number */
static void sort_seqctg_new(SeqCtg** seq_ctgs, int size)
{
  qsort(seq_ctgs, size, sizeof(SeqCtg *),
    (int (*)(const void*, const void*))seqctg_compare_new);
}

static gint
seqctg_compare_median(const SeqCtg** c1, const SeqCtg** c2) 
{

  if((*c1)->median < (*c2)->median) 
    return -1;
  else if((*c1)->median > (*c2)->median) 
    return 1;
  else
    return 0;
}

static gint
seqctg_compare_new(const SeqCtg** c1, const SeqCtg** c2) 
{
  if((*c1)->new_seqctg < (*c2)->new_seqctg) 
    return -1;
  else if((*c1)->new_seqctg > (*c2)->new_seqctg) 
    return 1;
  else
    return 0;
}

/*                   DEF: sort_coord 
 *Sort clones in ascending order of relevent coordinates*/
static void sort_coord(int* coord, int size)
{
  qsort(coord, size, sizeof(int),
    (int (*)(const void*, const void*))coord_compare_asce);
}

static gint
coord_compare_asce(const int* c1, const int* c2) 
{

  if(*c1 < *c2) 
    return -1;
  else if(*c1 > *c2) 
    return 1;
  else
    return 0;
}
